---
title: "Public Perceptions of Corporate Position-Taking on Abortion and Transgender Rights"
author: "Wayde Z.C. Marsh and Jordan Carr Peterson"
date: "12/19/22"
output:
  html_document:
    toc: true
    code_folding: hide
  latex_engine: xelatex
  number_sections: yes
mainfont: Helvetica
editor_options: 
  chunk_output_type: inline
---

# Loading packages and data

```{r, message = FALSE, warning = FALSE, results ='hide'}

##### Packages #####

library(descr)
library(ggplot2)
library(ggpubr)
library(foreign)
library(Rmisc)
library(plyr)
library(ggpubr)
library(ggeffects)
library(dplyr)
library(margins)
library(modelsummary)

knitr::opts_knit$set(root.dir = '~/Dropbox/Projects/CPT/Data')

```

## Load data

```{r}

setwd('~/Dropbox/Projects/CPT/Data') # set working directory

load('dat.Rda')

dat_1 <- subset(dat, study == 1)
dat_2 <- subset(dat, study == 2)
dat_3 <- subset(dat, study == 3)

```

# Modeling effects

## Policy preference DVs

### Anti-Trans Policies

```{r}

## Study I

m1_1r <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 1 & exp == 1 & gop == 1))
summary(m1_1r)

m1_1d <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 1 & exp == 1 & gop == 0))
summary(m1_1d)

## Study II

m1_2r <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 2 & exp == 1 & gop == 1))
summary(m1_2r)

m1_2d <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 2 & exp == 1 & gop == 0))
summary(m1_2d)

## Study III

m1_3r <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 3 & exp == 1 & gop == 1))
summary(m1_3r)

m1_3d <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 3 & exp == 1 & gop == 0))
summary(m1_3d)

```

### Anti-Abortion Policies

```{r}

## Study I

m1r_1r <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 1 & exp == 2 & gop == 1))
summary(m1r_1r)

m1r_1d <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 1 & exp == 2 & gop == 0))
summary(m1r_1d)

## Study II

m1r_2r <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 2 & exp == 2 & gop == 1))
summary(m1r_2r)

m1r_2d <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 2 & exp == 2 & gop == 0))
summary(m1r_2d)

## Study III

m1r_3r <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 3 & exp == 2 & gop == 1))
summary(m1r_3r)

m1r_3d <- lm(dv1 ~ factor(tr), dat = subset(dat, study == 3 & exp == 2 & gop == 0))
summary(m1r_3d)

```

### Tables

```{r}

modelsummary(models = list('Study I'  = m1_1r,
                           'Study II' = m1_2r,
                           'Study III' = m1_3r,
                           'Study I'  = m1_1d,
                           'Study II' = m1_2d,
                           'Study III' = m1_3d),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Study I'  = m1r_1r,
                           'Study II' = m1r_2r,
                           'Study III' = m1r_3r,
                           'Study I'  = m1r_1d,
                           'Study II' = m1r_2d,
                           'Study III' = m1r_3d),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

```

```{r}

# policy preferences

coef1 <- c(m1_1r$coefficients[2], m1_2r$coefficients[2], # reps coefs trans
           m1_3r$coefficients[2],
           m1_1d$coefficients[2], m1_2d$coefficients[2], # dem coefs trans
           m1_3d$coefficients[2],
           m1r_1r$coefficients[2], m1r_2r$coefficients[2], # reps coefs abort
           m1r_3r$coefficients[2],
           m1r_1d$coefficients[2], m1r_3d$coefficients[2], # dem coefs abort
           m1r_3d$coefficients[2]) 
se1 <- sqrt(diag(vcov(m1_1r))); se2 <- sqrt(diag(vcov(m1_2r))); se3 <- sqrt(diag(vcov(m1_3r))); # reps SEs trans
se4 <- sqrt(diag(vcov(m1_1d))); se5 <- sqrt(diag(vcov(m1_2d))); se6 <- sqrt(diag(vcov(m1_3d)));  # dems SE trans
se7 <- sqrt(diag(vcov(m1r_1r))); se8 <- sqrt(diag(vcov(m1r_2r))); se9 <- sqrt(diag(vcov(m1r_3r))); # reps SEs abort
se10 <- sqrt(diag(vcov(m1r_1d))); se11 <- sqrt(diag(vcov(m1r_2d))); se12 <- sqrt(diag(vcov(m1r_3d)));  # dems SE abort

hi1 <- c(coef1[1]+(1.96*se1[2]),   coef1[2]+(1.96*se2[2]),
         coef1[3]+(1.96*se3[2]),   coef1[4]+(1.96*se4[2]),
         coef1[5]+(1.96*se5[2]),   coef1[6]+(1.96*se6[2]),
         coef1[7]+(1.96*se7[2]),   coef1[8]+(1.96*se8[2]),
         coef1[9]+(1.96*se9[2]),   coef1[10]+(1.96*se10[2]),
         coef1[11]+(1.96*se11[2]), coef1[12]+(1.96*se12[2]))

lo1 <- c(coef1[1]-(1.96*se1[2]),   coef1[2]-(1.96*se2[2]),
         coef1[3]-(1.96*se3[2]),   coef1[4]-(1.96*se4[2]),
         coef1[5]-(1.96*se5[2]),   coef1[6]-(1.96*se6[2]),
         coef1[7]-(1.96*se7[2]),   coef1[8]-(1.96*se8[2]),
         coef1[9]-(1.96*se9[2]),   coef1[10]-(1.96*se10[2]),
         coef1[11]-(1.96*se11[2]), coef1[12]-(1.96*se12[2]))

party <- c('Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic')
policy <- c('Transgender Sports', 'Transgender Sports', 'Transgender Sports',
            'Transgender Sports', 'Transgender Sports', 'Transgender Sports',
            'Abortion', 'Abortion', 'Abortion',
            'Abortion', 'Abortion', 'Abortion')

study <- rep(c('Study I', 'Study II', 'Study III'), 4)

df_1 <- as.data.frame(cbind(coef1, hi1, lo1, party, policy, study))
df_1[, 1:3] <- lapply(df_1[, 1:3], as.numeric)
df_1[, 4:6] <- lapply(df_1[, 4:6], as.factor)

plot1 <- ggplot(df_1, aes(shape = study)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = policy, xmin = lo1,
                         xmax = hi1),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-20, 50) +
  scale_shape_manual(name = 'Study',
                    labels = c('Study I', 'Study II', 'Study III'),
                    values = c(1, 17, 12),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Policy Treatment',
                   labels = c('Transgender Sports',
                              'Abortion')) +
  scale_x_continuous("Estimated Effect on Policy Support \n(High Values Indicate more ``Liberal'' Position)") +
  labs(title = 'Corporate Position-Taking does not Impact Policy Preferences') +
  facet_wrap(~party) +
  # caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nError bars represent 95/% confidence intervals. Full results table in Supplementary Information.' +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot1

# ggsave('plot1.png', plot1, 'png', '~/Dropbox/Projects/CPT',
#     width = 190, height = 127, units = 'mm')

```

## CPT DVs

```{r}

### Support for Corporate Engagement in Politics

m2_1 <- lm(dv2 ~ factor(tr)*gop, dat = subset(dat, study == 1 & exp == 1))
summary(m2_1)

m2_2 <- lm(dv2 ~ factor(tr)*gop, dat = subset(dat, study == 2 & exp == 1))
summary(m2_2)

m2_3 <- lm(dv2 ~ factor(tr)*gop, dat = subset(dat, study == 3 & exp == 1))
summary(m2_3)

#### Breaking down by PID

m2_1r <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 1 & exp == 1 & gop == 1))
summary(m2_1r)

m2_1d <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 1 & exp == 1 & gop == 0))
summary(m2_1d)

m2_2r <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 2 & exp == 1 & gop == 1))
summary(m2_2r)

m2_2d <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 2 & exp == 1 & gop == 0))
summary(m2_2d)

m2_3r <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 3 & exp == 1 & gop == 1))
summary(m2_3r)

m2_3d <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 3 & exp == 1 & gop == 0))
summary(m2_3d)

### Anti-Abortion Policies

m2r_1 <- lm(dv2 ~ factor(tr)*gop, dat = subset(dat, study == 1 & exp == 2))
summary(m2r_1)

m2r_2 <- lm(dv2 ~ factor(tr)*gop, dat = subset(dat, study == 2 & exp == 2))
summary(m2r_2)

m2r_3 <- lm(dv2 ~ factor(tr)*gop, dat = subset(dat, study == 3 & exp == 2))
summary(m2r_3)

#### Breaking down by PID

m2r_1r <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 1 & exp == 2 & gop == 1))
summary(m2r_1r)

m2r_1d <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 1 & exp == 2 & gop == 0))
summary(m2r_1d)

m2r_2r <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 2 & exp == 2 & gop == 1))
summary(m2r_2r)

m2r_2d <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 2 & exp == 2 & gop == 0))
summary(m2r_2d)

m2r_3r <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 3 & exp == 2 & gop == 1))
summary(m2r_3r)

m2r_3d <- lm(dv2 ~ factor(tr), dat = subset(dat, study == 3 & exp == 2 & gop == 0))
summary(m2r_3d)

```

### Tables

```{r}

modelsummary(models = list('Study I'  = m2_1r,
                           'Study II' = m2_2r,
                           'Study III' = m2_3r,
                           'Study I'  = m2_1d,
                           'Study II' = m2_2d,
                           'Study III' = m2_3d),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Study I'  = m2r_1r,
                           'Study II' = m2r_2r,
                           'Study III' = m2r_3r,
                           'Study I'  = m2r_1d,
                           'Study II' = m2r_2d,
                           'Study III' = m2r_3d),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

```

```{r}

# CPT public opinion

coef2 <- c(m2_1r$coefficients[2], m2_2r$coefficients[2], # reps coefs trans
           m2_3r$coefficients[2],
           m2_1d$coefficients[2], m2_2d$coefficients[2], # dem coefs trans
           m2_3d$coefficients[2],
           m2r_1r$coefficients[2], m2r_2r$coefficients[2], # reps coefs abort
           m2r_3r$coefficients[2],
           m2r_1d$coefficients[2], m2r_3d$coefficients[2], # dem coefs abort
           m2r_3d$coefficients[2]) 
se1 <- sqrt(diag(vcov(m2_1r))); se2 <- sqrt(diag(vcov(m2_2r))); se3 <- sqrt(diag(vcov(m2_3r))); # reps SEs trans
se4 <- sqrt(diag(vcov(m2_1d))); se5 <- sqrt(diag(vcov(m2_2d))); se6 <- sqrt(diag(vcov(m2_3d)));  # dems SE trans
se7 <- sqrt(diag(vcov(m2r_1r))); se8 <- sqrt(diag(vcov(m2r_2r))); se9 <- sqrt(diag(vcov(m2r_3r))); # reps SEs abort
se10 <- sqrt(diag(vcov(m2r_1d))); se11 <- sqrt(diag(vcov(m2r_2d))); se12 <- sqrt(diag(vcov(m2r_3d)));  # dems SE abort

hi2 <- c(coef2[1]+(1.96*se1[2]),   coef2[2]+(1.96*se2[2]),
         coef2[3]+(1.96*se3[2]),   coef2[4]+(1.96*se4[2]),
         coef2[5]+(1.96*se5[2]),   coef2[6]+(1.96*se6[2]),
         coef2[7]+(1.96*se7[2]),   coef2[8]+(1.96*se8[2]),
         coef2[9]+(1.96*se9[2]),   coef2[10]+(1.96*se10[2]),
         coef2[11]+(1.96*se11[2]), coef2[12]+(1.96*se12[2]))

lo2 <- c(coef2[1]-(1.96*se1[2]),   coef2[2]-(1.96*se2[2]),
         coef2[3]-(1.96*se3[2]),   coef2[4]-(1.96*se4[2]),
         coef2[5]-(1.96*se5[2]),   coef2[6]-(1.96*se6[2]),
         coef2[7]-(1.96*se7[2]),   coef2[8]-(1.96*se8[2]),
         coef2[9]-(1.96*se9[2]),   coef2[10]-(1.96*se10[2]),
         coef2[11]-(1.96*se11[2]), coef2[12]-(1.96*se12[2]))

party <- c('Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic')
policy <- c('Transgender Sports', 'Transgender Sports', 'Transgender Sports',
            'Transgender Sports', 'Transgender Sports', 'Transgender Sports',
            'Abortion', 'Abortion', 'Abortion',
            'Abortion', 'Abortion', 'Abortion')

study <- rep(c('Study I', 'Study II', 'Study III'), 4)

df_2 <- as.data.frame(cbind(coef2, hi2, lo2, party, policy, study))
df_2[, 1:3] <- lapply(df_2[, 1:3], as.numeric)
df_2[, 4:6] <- lapply(df_2[, 4:6], as.factor)

plot2 <- ggplot(df_2, aes(shape = study)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef2, y = policy, xmin = lo2,
                         xmax = hi2),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-20, 50) +
  scale_shape_manual(name = 'Study',
                    labels = c('Study I', 'Study II', 'Study III'),
                    values = c(1, 17, 12),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Policy Treatment',
                   labels = c('Transgender Sports',
                              'Abortion')) +
  scale_x_continuous("Estimated Effect on Support for Corporate Political Engagement") +
  labs(title = 'Corporate Position-Taking Impacts Public Opinion \non Corporate Political Engagement ') +
  facet_wrap(~party) +
  # caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nError bars represent 95/% confidence intervals. Full results table in Supplementary Information.' +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot2

# ggsave('plot2.png', plot2, 'png', '~/Dropbox/Projects/CPT',
#     width = 190, height = 127, units = 'mm')

```

## Secondary Experiment

## Secondary Experiment--Punishment and Reward

```{r}

## Punish v. Reward CPT

m3_2 <- lm(dv3_reps ~ tr, dat = subset(dat, study == 2))
summary(m3_2)

m3_3 <- lm(dv3_reps ~ tr, dat = subset(dat, study == 3))
summary(m3_3)

#### Breaking down by PID

m3_1r <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 1 & exp == 1 & gop == 1))
summary(m3_1r)

m3_1d <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 1 & exp == 1 & gop == 0))
summary(m3_1d)

m3_2r <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 2 & exp == 1 & gop == 1))
summary(m3_2r)

m3_2d <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 2 & exp == 1 & gop == 0))
summary(m3_2d)

m3_3r <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 3 & exp == 1 & gop == 1))
summary(m3_3r)

m3_3d <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 3 & exp == 1 & gop == 0))
summary(m3_3d)


m3r_1r <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 1 & exp == 2 & gop == 1))
summary(m3r_1r)

m3r_1d <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 1 & exp == 2 & gop == 0))
summary(m3r_1d)

m3r_2r <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 2 & exp == 2 & gop == 1))
summary(m3r_2r)

m3r_2d <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 2 & exp == 2 & gop == 0))
summary(m3r_2d)

m3r_3r <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 3 & exp == 2 & gop == 1))
summary(m3r_3r)

m3r_3d <- lm(dv3 ~ factor(exp2), dat = subset(dat, study == 3 & exp == 2 & gop == 0))
summary(m3r_3d)

## Effect of abortion v. transgender sports

m4_1r <- lm(dv3 ~ factor(exp), dat = subset(dat, study == 1 & gop == 1))
summary(m4_1r)

m4_1d <- lm(dv3 ~ factor(exp), dat = subset(dat, study == 1 & gop == 0))
summary(m4_1d)

m4_2r <- lm(dv3 ~ factor(exp), dat = subset(dat, study == 2 & gop == 1))
summary(m4_2r)

m4_2d <- lm(dv3 ~ factor(exp), dat = subset(dat, study == 2 & gop == 0))
summary(m4_2d)

m4_3r <- lm(dv3 ~ factor(exp), dat = subset(dat, study == 3 & gop == 1))
summary(m4_3r)

m4_3d <- lm(dv3 ~ factor(exp), dat = subset(dat, study == 3 & gop == 0))
summary(m4_3d)

## Effect of CPT treatment

m3_2.2r <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 2 & exp == 1 & gop == 1))
summary(m3_2.2r)

m3_2.2d <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 2 & exp == 1 & gop == 0))
summary(m3_2.2d)

m3_3.2r <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 3 & exp == 1 & gop == 1))
summary(m3_3.2r)

m3_3.2d <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 3 & exp == 1 & gop == 0))
summary(m3_3.2d)

m3_2r.2r <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 2 & exp == 2 & gop == 1))
summary(m3_2r.2r)

m3_2r.2d <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 2 & exp == 2 & gop == 0))
summary(m3_2r.2d)

m3_3r.2r <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 3 & exp == 2 & gop == 1))
summary(m3_3r.2r)

m3_3r.2d <- lm(dv3 ~ factor(exp2)*factor(tr), dat = subset(dat, study == 3 & exp == 2 & gop == 0))
summary(m3_3r.2d)

```

### Tables

```{r}

modelsummary(models = list('Study I'  = m3_1r,
                           'Study II' = m3_2r,
                           'Study III' = m3_3r,
                           'Study I'  = m3_1d,
                           'Study II' = m3_2d,
                           'Study III' = m3_3d),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Study I'  = m3r_1r,
                           'Study II' = m3r_2r,
                           'Study III' = m3r_3r,
                           'Study I'  = m3r_1d,
                           'Study II' = m3r_2d,
                           'Study III' = m3r_3d),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

```

```{r}

# CPT punishment 

coef3.1 <- c(m3_1r$coefficients[2], m3_2r$coefficients[2], # reps coefs trans
           m3_3r$coefficients[2],
           m3_1d$coefficients[2], m3_2d$coefficients[2], # dem coefs trans
           m3_3d$coefficients[2],
           m3r_1r$coefficients[2], m3r_2r$coefficients[2], # reps coefs abort
           m3r_3r$coefficients[2],
           m3r_1d$coefficients[2], m3r_3d$coefficients[2], # dem coefs abort
           m3r_3d$coefficients[2]) 
se1 <- sqrt(diag(vcov(m3_1r))); se2 <- sqrt(diag(vcov(m3_2r))); se3 <- sqrt(diag(vcov(m3_3r))); # reps SEs trans
se4 <- sqrt(diag(vcov(m3_1d))); se5 <- sqrt(diag(vcov(m3_2d))); se6 <- sqrt(diag(vcov(m3_3d)));  # dems SE trans
se7 <- sqrt(diag(vcov(m3r_1r))); se8 <- sqrt(diag(vcov(m3r_2r))); se9 <- sqrt(diag(vcov(m3r_3r))); # reps SEs abort
se10 <- sqrt(diag(vcov(m3r_1d))); se11 <- sqrt(diag(vcov(m3r_2d))); se12 <- sqrt(diag(vcov(m3r_3d)));  # dems SE abort

hi2 <- c(coef3.1[1]+(1.96*se1[2]),   coef3.1[2]+(1.96*se2[2]),
         coef3.1[3]+(1.96*se3[2]),   coef3.1[4]+(1.96*se4[2]),
         coef3.1[5]+(1.96*se5[2]),   coef3.1[6]+(1.96*se6[2]),
         coef3.1[7]+(1.96*se7[2]),   coef3.1[8]+(1.96*se8[2]),
         coef3.1[9]+(1.96*se9[2]),   coef3.1[10]+(1.96*se10[2]),
         coef3.1[11]+(1.96*se11[2]), coef3.1[12]+(1.96*se12[2]))

lo2 <- c(coef3.1[1]-(1.96*se1[2]),   coef3.1[2]-(1.96*se2[2]),
         coef3.1[3]-(1.96*se3[2]),   coef3.1[4]-(1.96*se4[2]),
         coef3.1[5]-(1.96*se5[2]),   coef3.1[6]-(1.96*se6[2]),
         coef3.1[7]-(1.96*se7[2]),   coef3.1[8]-(1.96*se8[2]),
         coef3.1[9]-(1.96*se9[2]),   coef3.1[10]-(1.96*se10[2]),
         coef3.1[11]-(1.96*se11[2]), coef3.1[12]-(1.96*se12[2]))

party <- c('Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic',
           'Republican', 'Republican', 'Republican',
           'Democratic', 'Democratic', 'Democratic')
policy <- c('Transgender Sports', 'Transgender Sports', 'Transgender Sports',
            'Transgender Sports', 'Transgender Sports', 'Transgender Sports',
            'Abortion', 'Abortion', 'Abortion',
            'Abortion', 'Abortion', 'Abortion')

study <- rep(c('Study I', 'Study II', 'Study III'), 4)

df_3.1 <- as.data.frame(cbind(coef3.1, hi2, lo2, party, policy, study))
df_3.1[, 1:3] <- lapply(df_3.1[, 1:3], as.numeric)
df_3.1[, 4:6] <- lapply(df_3.1[, 4:6], as.factor)

plot3.1 <- ggplot(df_3.1, aes(shape = study)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef3.1, y = policy, xmin = lo2,
                         xmax = hi2),
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  xlim(-20, 50) +
  scale_shape_manual(name = 'Study',
                    labels = c('Study I', 'Study II', 'Study III'),
                    values = c(1, 17, 12),
                    guide = guide_legend(reverse = T)) +
  scale_y_discrete(name = 'Policy Treatment',
                   labels = c('Transgender Sports',
                              'Abortion')) +
  scale_x_continuous("Estimated Effect of Legislator Action on Support for Corporate Postion-Taking Constraint") +
  labs(title = 'Information on Partisan Reward/Punishment by Legislators Impacts Public Opinion \non Corporate Position-Taking Constraints ') +
  facet_wrap(~party) +
  # caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identification). \nError bars represent 95/% confidence intervals. Full results table in Supplementary Information.' +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot3.1

# ggsave('plot3.1.png', plot3.1, 'png', '~/Dropbox/Projects/CPT',
#     width = 190, height = 127, units = 'mm')


```


